%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% O-Ring Production Networks
%   Author: Demir, Fieler, Xu, Yang
%   Last updated: January 2021
%
% Purpose: Standard Error (Full Model)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all; diary off; clc; clear; 


%% Parameters

    % load point estimates
    load('../../output/estimates/fullmodel_estimates.mat', 'est');            
    theta = est;

    % parameters
    define_param
    
    % increase upper bound
    param.ub = 10;
    
    
%% Jacobian (G)

    % initialize
    Nmom = 38;
    Npar = size(theta,2);
    G = zeros(Nmom,Npar);
    iter_inner_h = zeros(1,Npar); iter_outer_h = zeros(1,Npar);
    iter_inner_l = zeros(1,Npar); iter_outer_l = zeros(1,Npar);
    
    % perturbation
    perturb = 0.01;
    h = abs(theta)*perturb;

    % jacobian
    for i = 1:Npar  
        
        fprintf('\nParameter %.0f\n', i);
        theta_h = theta; theta_h(i) = theta(i)+h(i);
        theta_l = theta; theta_l(i) = theta(i)-h(i);
        [moments_h, iter_h] = sim_moments(theta_h, param);
        [moments_l, iter_l] = sim_moments(theta_l, param);
        G(:,i) = (moments_h-moments_l)'./(2*h(i));        
        iter_inner_h(i) = iter_h(1); iter_outer_h(i) = iter_h(2);
        iter_inner_l(i) = iter_l(1); iter_outer_l(i) = iter_l(2);
        
    end    
    
    % save
    mindex = [repelem([8:-1:3],5), repelem([2:-1:1],4)]';
    quintile = [repmat([1:5],1,6), repmat([2:5],1,2)]';
    Jacobian = [mindex, quintile, G];
    mkdir('../../output/temp') 
    writematrix(Jacobian, '../../output/temp/jacobian.csv');    
    

 
